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Extinction  appears  ubiquitously  in  many  fields,  including  chemical  reactions,  population 
biology,  evolution  and  epidemiology.  Even  though  extinction  as  a  random  process  is  a  rare 
event,  its  occurrence  is  observed  in  large  finite  populations.  Extinction  occurs  when  fluctu¬ 
ations  owing  to  random  transitions  act  as  an  effective  force  that  drives  one  or  more 
components  or  species  to  vanish.  Although  there  are  many  random  paths  to  an  extinct 
state,  there  is  an  optimal  path  that  maximizes  the  probability  to  extinction.  In  this  paper, 
we  show  that  the  optimal  path  is  associated  with  the  dynamical  systems  idea  of  having  maxi¬ 
mum  sensitive  dependence  to  initial  conditions.  Using  the  equivalence  between  the  sensitive 
dependence  and  the  path  to  extinction,  we  show  that  the  dynamical  systems  picture  of  extinc¬ 
tion  evolves  naturally  towards  the  optimal  path  in  several  stochastic  models  of  epidemics. 
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1.  INTRODUCTION 

Determining  the  conditions  for  epidemic  extinction  is 
an  important  public  health  problem.  Global  eradication 
of  an  infectious  disease  has  rarely  been  achieved,  but  it 
continues  to  be  a  public  health  goal  for  polio  [1]  and 
many  other  diseases,  including  childhood  diseases. 
More  commonly,  disease  extinction,  or  fade  out,  is 
local  and  may  be  followed  by  a  reintroduction  of  the  dis¬ 
ease  from  other  regions  [2,3].  Extinction  may  also  occur 
for  individual  strains  of  a  multi-strain  disease  [4],  such 
as  influenza  or  dengue  fever.  Since  extinction  occurs 
in  finite  populations,  it  depends  critically  on  local  com¬ 
munity  size  [5-8].  Moreover,  it  is  important  to  know 
how  parameters  affect  the  chance  of  extinction  for  pre¬ 
dicting  the  dynamics  of  outbreaks  and  for  developing 
control  strategies  to  promote  epidemic  extinction  [9]. 
The  determination  of  extinction  risk  is  also  of  interest 
in  related  fields,  such  as  evolution  and  ecology.  For 
example,  in  the  neutral  theory  of  ecology,  bio-diversity 
arises  from  the  interplay  between  the  introduction  and 
extinction  of  species  [10,11]. 

In  general,  extinction  occurs  in  discrete,  finite  popu¬ 
lations  undergoing  stochastic  effects  owing  to  random 
transitions  or  perturbations.  The  origins  of  stochasticity 
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may  be  internal  to  the  system  or  may  arise  from  the  exter¬ 
nal  environment.  Small  population  size,  low  contact 
frequency  for  frequency-dependent  transmission,  compe¬ 
tition  for  resources  and  evolutionary  pressure  [12],  as 
well  as  heterogeneity  in  populations  and  transmission 
[13] ,  may  all  be  determining  factors  for  extinction  to  occur. 

Extinction  risk  is  affected  by  the  nature  and  strength 
of  the  noise  [14],  as  well  as  other  factors,  including  out¬ 
break  amplitude  [15]  and  seasonal  phase  occurrence 
[16].  For  large  populations,  the  intensity  of  internal 
population  noise  is  generally  small.  However,  a  rare, 
large  fluctuation  can  occur  with  non-zero  probability 
and  the  system  may  be  able  to  reach  the  extinct  state. 
For  respiratory  diseases  such  as  severe  acute  respiratory 
syndrome  (SARS),  super-spreading  may  account  for 
these  large  stochastic  fluctuations  [17].  Since  the  extinct 
state  is  absorbing  owing  to  effective  stochastic  forces, 
eventual  extinction  is  guaranteed  when  there  is  no 
source  of  reintroduction  [18-20].  However,  because 
fade  outs  are  usually  rare  events  in  large  populations, 
typical  timescales  for  extinction  may  be  extremely  long. 

Stochastic  population  models  of  finite  populations, 
which  include  extinction  processes,  are  effectively 
described  using  the  master  equation  formalism.  Sto¬ 
chastic  master  equations  are  commonly  used  in 
statistical  physics  when  dealing  with  chemical  reaction 
processes  [21]  and  predict  probabilities  of  rare  events 
occurring  at  certain  times.  For  many  problems  invol¬ 
ving  extinction  in  large  populations,  if  the  probability 
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distribution  of  the  population  is  approximately  station¬ 
ary,  the  probability  of  extinction  is  a  function  that 
decreases  exponentially  with  increasing  population 
size.  The  exponent  in  this  function  scales  as  a  determi¬ 
nistic  quantity  called  the  action  [22].  It  can  be  shown 
that  a  trajectory  that  brings  the  system  to  extinction 
is  very  likely  to  lie  along  a  most  probable  path  called 
the  optimal  extinction  trajectory  or  optimal  path.  It  is 
a  remarkable  property  that  a  deterministic  quantity 
such  as  the  action  can  predict  the  probability  of  extinc¬ 
tion,  which  is  inherently  a  stochastic  process  [9,23]. 

Locating  the  optimal  path  is  desirable  because  the 
quantity  of  interest,  the  extinction  rate,  depends  on 
the  probability  to  traverse  this  path,  and  the  effect  of 
a  control  strategy  on  extinction  rate  can  be  determined 
by  its  effect  on  the  optimal  path  [9].  By  employing  an 
optimal  path  formalism,  we  convert  the  stochastic  pro¬ 
blem  to  a  mechanistic  dynamical  systems  problem. 
In  contrast  to  approaches  based  on  diffusive  processes, 
which  are  valid  only  in  the  limit  of  large  system  sizes 
[24,25],  this  dynamical  systems  approach  can  give  accu¬ 
rate  estimates  for  the  extinction  time  even  for  small 
populations  if  the  action  is  sufficiently  large.  Addition¬ 
ally,  unlike  other  methods  that  are  used  to  estimate 
lifetimes,  this  approach  enables  one  both  to  estimate  life¬ 
times  and  to  draw  conclusions  about  the  path  taken  to 
extinction.  This  more  detailed  understanding  of  how 
extinction  occurs  may  lead  to  new  stochastic  control 
strategies  [9]. 

In  this  paper,  we  show  that  locating  the  optimal 
extinction  trajectory  can  be  achieved  naturally  by  evol¬ 
ving  a  dynamical  system  that  converges  to  the  optimal 
path.  The  method  is  based  on  computing  finite-time 
Lyapunov  exponents  (FTLE),  which  have  previously 
been  used  to  find  coherent  structures  in  fluid  flows 
[26-31].  The  FTLE  provides  a  measure  of  how  sensi¬ 
tively  the  system’s  future  behaviour  depends  on  its 
current  state.  We  argue  that  the  system  displays  maxi¬ 
mum  sensitivity  near  the  optimal  extinction  trajectory, 
which  enables  us  to  dynamically  evolve  towards  the 
optimal  escape  trajectory  using  FTLE  calculations. 
For  several  models  of  epidemics  that  contain  internal 
or  external  noise  sources,  we  illustrate  the  power  of 
our  method  to  locate  optimal  extinction  trajectories. 
Although  our  examples  are  taken  from  infectious  dis¬ 
ease  models,  the  approach  is  general  and  is  applicable 
to  any  extinction  process  or  escape  process. 


2.  STOCHASTIC  MODELLING  IN  THE 
LARGE  POPULATION  LIMIT 

To  introduce  our  idea  of  dynamically  constructing  an  opti¬ 
mal  path  to  extinction  in  stochastic  systems,  we  show 
its  application  to  a  stochastic  susceptible- infectious  - 
recovered  (SIR)  epidemiological  model.  Details  of  the 
SIR  model  can  be  found  in  the  electronic  supplementary 
material,  §1.  Figure  1  shows  the  probability  density  of 
extinction  prehistory  in  the  susceptible -infectious  (SI) 
plane.  The  probability  density  was  numerically  computed 
using  20  000  stochastic  SIR  trajectories  that  ended  in 
extinction.  Trajectories  are  aligned  at  their  extinction 
point.  From  the  extinction  point,  the  prehistory  of  each 
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susceptibles 

Figure  1.  Probability  density  of  extinction  prehistory  and  the 
optimal  path  to  extinction  for  the  stochastic  SIR  epidemiological 
model.  Colours  indicate  the  probability  density  (the  colour  bar 
values  have  been  normalized,  with  lighter  colours  corresponding 
to  higher  probability)  for  20  000  stochastic  realizations.  The 
results  were  computed  using  Monte  Carlo  simulations,  and  details 
of  the  sampling  of  the  trajectories  are  described  in  the  text.  The 
curve  is  the  numerically  predicted  optimal  path  to  extinction. 
Note  that  the  optimal  path  to  extinction  lies  on  the  peak  of  the 
probability  density  of  extinction  prehistory.  The  population  is 
3  x  106  individuals,  with  R0  «  15  (contact  rate  (3  =  1500, 
recovery  rate  y  =  100  and  birth- death  rate  p  =  0.2). 


trajectory  up  to  the  last  outbreak  of  infection  is  con¬ 
sidered.  Small  fluctuations  in  the  infectious  population 
are  not  considered  in  identifying  the  last  outbreak.  In 
this  way,  we  restrict  the  analysis  to  the  interval  between 
the  last  large  outbreak  of  infection  and  the  extinction 
point.  The  resulting  (S,I)  pairs  of  susceptible  and  infec¬ 
tious  individuals  are  then  binned  and  plotted  in  the  SI 
plane  [32] .  The  resulting  discrete  density  has  been  colour 
coded  so  that  the  brighter  regions  correspond  to  higher 
density  of  trajectories.  The  figure  shows  that,  among  all 
the  paths  that  the  stochastic  system  can  take  to  reach 
the  extinct  state,  there  is  one  path  that  has  the  highest 
probability  of  occurring.  This  is  the  optimal  path  to 
extinction.  One  can  see  that  the  optimal  path  to  extinction 
lies  on  the  peak  of  the  probability  density  of  the  extinction 
prehistory.  It  should  be  noted  that  extinction  for  the 
stochastic  SIR  model  has  been  studied  previously  [33] . 

The  optimal  path  can  be  obtained  using  methods  of 
statistical  physics.  In  figure  1,  the  numerical  prediction 
of  the  entire  optimal  trajectory  for  the  stochastic  SIR 
system  has  been  overlaid  on  the  probability  density  of 
extinction  prehistory  that  was  found  using  stochastic 
simulation.  The  trajectory  spirals  away  from  the  ende¬ 
mic  state,  with  larger  and  larger  oscillations  until  it 
hits  the  extinct  state.  The  agreement  between  the  sto¬ 
chastically  simulated  optimal  path  to  extinction  and 
the  predicted  optimal  path  is  excellent. 

The  curve  of  figure  1  is  obtained  by  recasting  the  sto¬ 
chastic  problem  in  a  deterministic  form.  The  evolution 
of  the  probability  of  finding  a  stochastic  system  in  a 
given  state  X  at  a  given  time  t  is  described  by  the 
master  equation  [34].  Solving  the  master  equation 
would  provide  a  complete  description  of  the  time  evol¬ 
ution  of  the  stochastic  system,  but  in  general  it  is  a 
difficult  task  to  obtain  explicit  solutions  for  the 
master  equation.  Thus,  one  generally  resorts  to  approxi¬ 
mations  to  the  solution;  i.e.  one  considers  an  ansatz  for 
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Figure  2.  (a)  Schematic  showing  the  fixed  points  and  heteroclinic  trajectories  (trajectories  connecting  fixed  points).  Coordinate 
axes  are  dashed  lines.  The  noise  coordinate  is  indicated  by  the  momentum  ( p)  coordinate.  The  deterministic  manifold  ( p  =  0)  is 
indicated  in  blue.  Deterministically,  the  extinct  state  is  repelling  and  endemic  state  is  attracting.  However,  the  endemic  state  has 
unstable  directions  for  non-zero  noise  (p  #  0),  and  the  optimal  path  (red)  is  the  trajectory  carrying  the  system  from  the  endemic 
state  to  a  stochastic  extinct  state.  ( b )  Schematic  showing  the  path  from  the  endemic  state  (blue)  to  the  extinct  state  (red).  The 
optimal  path  leaves  the  endemic  point  along  an  unstable  manifold,  and  connects  to  the  extinct  state  along  its  stable  manifold. 
The  magenta  and  green  dashed  lines  represent  trajectories  initially  separated  by  the  optimal  path.  The  initial  starting  distance 
between  trajectories  near  the  endemic  state  expands  exponentially  in  forward  time  (shown  by  the  cyan  lines).  Locally,  this  shows 
that  the  finite-time  Lyapunov  measure  of  sensitivity  with  respect  to  initial  data  is  maximal  along  the  optimal  path. 


the  probability  density.  In  this  case,  since  extinction 
of  finite  populations  is  a  rare  event,  we  will  be  interested 
in  examining  events  that  occur  in  the  tail  of  the 
probability  distribution.  Therefore,  the  distribution  is 
assumed  to  take  the  form: 

p(X,i)»exp(-JV5(*)),  (2.1) 

where  p(X,  t )  is  the  probability  density  of  finding  the 
system  in  the  state  X  at  time  t,  N  is  the  size  of  the 
population,  x  =  X/Ais  the  normalized  state  (e.g.  in  an 
epidemic  model,  the  fraction  of  the  population  in  the  var¬ 
ious  compartments),  and  S  is  a  deterministic  state 
function  known  in  classical  physics  as  the  action.  Equation 
(2.1)  describes  the  relationship  between  the  action  and  the 
probability  density  and  is  based  on  an  assumption  for 
how  the  probability  scales  with  the  population  size.  The 
action  is  the  negative  of  the  natural  log  of  the  stationary 
probability  distribution  divided  by  the  population  size. 
Therefore,  the  probability  (if  we  normalize  the  popu¬ 
lation)  is  roughly  given  by  the  exponential  of  the  action. 
Intuitively,  equation  (2.1)  expresses  the  assumption  that 
the  probability  of  occurrence  of  extreme  events,  such  as 
extinction,  lies  in  the  tails  of  the  probability  distribution, 
which  falls  steeply  away  from  the  steady  state. 

This  approximation  leads  to  a  conserved  quan¬ 
tity  that  is  called  the  Hamiltonian  [35].  From  the 
Hamiltonian,  one  can  find  a  set  of  conservative  ordinary 
differential  equations  (ODEs)  that  are  known  as 
Hamilton’s  equations.  These  ODEs  describe  the  time 
evolution  of  the  system  in  terms  of  state  variables  x, 
which  are  associated  with  the  population  in  each  com¬ 
partment.  For  the  SIR  example,  x  is  the  vector 
(S',/,  A).  In  addition  to  the  state  variables,  the  equations 
contain  conjugate  momenta  variables,  The  conju¬ 
gate  momenta,  or  noise,  account  for  the  uncertainty 
associated  with  the  system  being  in  a  given  state  at  a 
given  time  owing  to  the  stochastic  interactions  among 
the  individuals  of  the  population.  These  ODEs  can  be 


constructed  from  information  in  the  master  equation 
about  the  possible  transitions  and  transition  rates  in 
the  system.  Details  can  be  found  in  appendix  A. 

Integration  of  the  ODEs  with  the  appropriate  bound¬ 
ary  conditions  will  then  give  the  optimal  evolution  of  the 
system  under  the  influence  of  the  noise.  Boundary 
conditions  are  chosen  to  be  fixed  points  of  the  system. 
A  typical  case  is  shown  schematically  in  figure  2  a.  Deter¬ 
ministically,  the  endemic  state  is  attracting  and  the 
extinct  state  repelling.  However,  introducing  stochasti- 
city  allows  the  system  to  leave  the  deterministic 
manifold  along  an  unstable  direction  of  the  endemic 
state,  corresponding  to  non-zero  noise.  Stochasticity 
leads  to  an  additional  extinct  state  that  arises  owing  to 
the  general  non- Gaussian  nature  of  the  noise.  For  the 
extinction  process  of  figure  1,  boundary  conditions 
were  the  system  leaving  the  endemic  steady  state  and 
asymptotically  approaching  the  stochastic  extinct  state. 

In  general,  the  optimal  extinction  path  is  an  unstable 
dynamical  object,  and  this  reflects  extinction  as  a  rare 
event.  This  has  led  many  authors  to  consider  how 
extinction  rates  scale  with  respect  to  a  parameter 
close  to  a  bifurcation  point  [9,33,36,37],  where  the 
dynamics  are  very  slow.  For  an  epidemic  model,  this 
means  that  the  reproductive  rate  R0  should  be  greater 
than  but  very  close  to  1.  However,  most  real  diseases 
have  Rq  larger  than  1.5,  which  translates  into  a  faster 
growth  rate  from  the  extinct  state.  In  general,  in  order 
to  obtain  analytical  scaling  results,  one  must  obtain  the 
ODEs  for  the  optimal  path  either  analytically  (using 
the  classical  theory  of  large  fluctuations  mentioned 
within  this  section  and  described  in  detail  in  appendix 
A)  or  numerically  (using  shooting  methods  for  boundary 
value  problems).  This  task  may  be  impossible  or  extre¬ 
mely  cumbersome,  especially  when  the  system  is  far 
from  the  bifurcation  point.  In  the  following  section,  we 
demonstrate  how  to  evolve  naturally  to  the  optimal 
path  to  extinction  using  a  dynamical  systems  approach. 
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3.  FINITE-TIME  LYAPUNOV  EXPONENTS 

We  consider  a  velocity  field  that  is  defined  over  a  finite 
time  interval  and  is  given  by  Hamilton’s  equations  of 
motion.  Such  a  velocity  field,  when  starting  from  an 
initial  condition,  has  a  unique  solution.  The  continuous 
dynamical  system  has  quantities,  known  as  Lyapunov 
exponents,  which  are  associated  with  the  trajectory  of 
the  system  in  an  infinite  time  limit,  and  which  measure 
the  average  growth  rates  of  the  linearized  dynamics 
about  the  trajectory.  To  find  the  FTLE,  one  computes 
the  Lyapunov  exponents  on  a  restricted  finite  time  inter¬ 
val.  For  each  initial  condition,  the  exponents  provide  a 
measure  of  their  sensitivity  to  small  perturbations. 
Therefore,  the  FTLE  is  a  measure  of  the  local  sensitivity 
to  initial  data.  Details  regarding  the  FTLE  can  be  found 
in  the  electronic  supplementary  material,  §2. 

The  FTLE  measurements  can  be  shown  to  exhibit 
‘ridges’  of  local  maxima.  The  ridges  of  the  field  indicate 
the  location  of  attracting  (backward  time  FTLE  field) 
and  repelling  (forward  time  FTLE  field)  structures.  In 
two-dimensional  space,  the  ridge  is  a  curve,  which 
locally  maximizes  the  FTLE  field  so  that  transverse  to 
the  ridge  one  finds  the  FTLE  to  be  a  local  maximum. 
What  is  remarkable  is  that  the  FTLE  ridges  correspond 
to  the  optimal  path  trajectories,  which  we  heuristically 
argue  in  the  electronic  supplementary  material,  §3.  The 
basic  idea  is  that  since  the  optimal  path  is  inherently 
unstable,  the  FTLE  shows  that,  locally,  the  path  is 
also  the  most  sensitive  to  initial  data.  Figure  2b  shows 
a  schematic  that  demonstrates  why  the  optimal  path 
has  a  local  maximum  to  sensitivity.  If  one  chooses  an 
initial  point  on  either  side  of  the  path  near  the  endemic 
state,  the  two  trajectories  will  separate  exponentially 
in  time.  This  is  due  to  the  fact  that  both  extinct  and 
endemic  states  are  unstable,  and  the  connecting  trajec¬ 
tory  defining  the  path  is  unstable  as  well.  Any  initial 
points  starting  near  the  optimal  path  will  leave  the 
neighbourhood  in  short  time. 


4.  EVOLVING  TOWARDS  THE  OPTIMAL 
PATH  USING  FINITE-TIME  LYAPUNOV 
EXPONENTS 

We  now  apply  our  theory  of  dynamical  sensitivity  to 
the  problem  of  locating  optimal  paths  to  extinction 
for  several  examples.  We  consider  the  case  of  internal 
fluctuations,  where  the  noise  is  not  known  a  priori ,  as 
well  as  the  case  of  external  noise.  In  each  case,  the  inter¬ 
action  of  the  noise  and  state  of  the  system  begins  by 
finding  the  equations  of  motion  that  describe  the 
unstable  flow.  These  equations  of  motion  are  then  used 
to  compute  the  ridges  corresponding  to  maximum 
FTLE,  which  in  turn  correspond  to  the  optimal 
extinction  paths  [38] . 


4-.1.  Extinction  in  a  branching- annihilation 
process 

For  an  example  of  a  system  with  internal  fluctuations, 
which  has  an  analytical  solution,  consider  extinction 


in  the  stochastic  branching -annihilation  process 

A^2A  and  2  A-^0,  (4.1) 

where  A  and  p  >  0  are  constant  reaction  rates  [39,40]. 
Equation  (4.1)  is  a  single  species  birth- death  process 
and  can  be  thought  of  as  a  simplified  form  of  the  Verhulst 
logistic  model  for  population  growth  [41] .  The  mean  field 
equation  for  the  average  number  of  individuals  n  in  the 
infinite  population  limit  is  given  by  h  =  An  —  pn 2.  The 
stochastic  process  given  by  equation  (4.1)  contains 
intrinsic  noise  that  arises  from  the  randomness  of 
the  reactions  and  the  fact  that  the  population  consists 
of  discrete  individuals.  This  intrinsic  noise  can  generate 
a  rare  sequence  of  events  that  causes  the  system  to 
evolve  to  the  extinct  state.  The  probability  Pn{t)  to 
observe,  at  time  t,  n  individuals  is  governed  by  the 
master  equation 

Pn  =*•  ~  i(n  +  2 )(n  +  l)Pn+2  -  n(n  -  1  )Pn] 

+  A[(n-l)P„-1-nP„].  (4.2) 


The  Hamiltonian  associated  with  this  system  is 


H(q,p)  =  (A(l  +  p)  ~^(2  +  p)q)qp, 

(4.3) 

where  q  is  a  conjugate  coordinate  related  to  n  through  a 
transformation  [40] ,  and  p  plays  the  role  of  the  momen¬ 
tum.  The  equations  of  motion  are  given  by 

dH 

q =  ~~q~  =  tfW1  +  2p)  _  M1  +  p)q\ 

(4.4a) 

and 

dH 

p=  --g^  =  p[p-(2  +  p)q-K1  +  p)}- 

(4.46) 

The  mean  field  is  retrieved  in  equation  (4.4a)  when  p  = 
0  (no  fluctuations  or  noise).  The  Hamiltonian  has  three 
zero-energy  curves.  The  first  is  the  mean-field  zero- 
energy  line  p—  0  (no  fluctuations),  which  contains 
two  unstable  points  hi  =  (p,q)  =  (0,A//i)  and  h0  = 
( p,q )  =  (0,0).  The  second  is  the  extinction  line  q=  0 
(trivial  solution),  which  contains  another  unstable 
point  h2  =  (p,q)  =  (  — 1,0).  The  third  zero-energy 
curve  is  non-trivial  and  is 


2A(1  +  p) 

p{2  +  p) 


(4.5) 


The  segment  of  the  curve  given  by  equation  (4.5), 
which  lies  between  —  1  <  p  <  0  corresponds  to  a  (het¬ 
eroclinic)  trajectory  that  exits,  at  t  =  —  oo?  the  point 
hi  along  its  unstable  manifold  and  enters,  at  t=  oo? 
the  point  h2  along  its  stable  manifold.  This  trajectory 
is  the  optimal  path  to  extinction  and  describes  the 
most  probable  sequence  of  events,  which  evolves  the 
system  from  a  quasi-stationary  state  to  extinction  [40]. 

To  show  that  the  FTLE  evolves  to  the  optimal  path, 
we  calculate  the  FTLE  field  using  the  system  of  Hamil¬ 
ton’s  equations  given  by  equations  (4.4a,  b).  Figure  3a 
shows  both  the  forward  and  backward  FTLE  plot  com¬ 
puted  for  the  finite  time  T—  6,  with  A  =  2.0  and  p  = 
0.5.  In  this  example,  as  well  as  the  following  two 
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Figure  3.  (a)  Forward  and  backward  FTLE  flow  fields  for  the  branching -annihilation  process  computed  using  the  Hamiltonian 
given  by  equation  (4.3)  with  A  =  2.0  and  pi  =  0.5.  The  integration  time  is  T  =  6  with  an  integration  step  size  of  t  =  0.1  and  a  grid 
resolution  of  0.01  in  both  q  and  p  (momentum).  The  three  zero-energy  curves  are  given  by  the  ridges  of  maximal  FTLE  and  are 
overlaid  with  the  analytical  solution  of  these  curves  given  by  p  =  0,  q  =  0  and  equation  (4.5).  (6)  Forward  and  backward  FTLE 
flow  fields  for  the  SIS  epidemic  model  with  external  fluctuations.  The  flow  fields  were  computed  using  the  Lagrangian  given  by 
equation  (4.9)  with  (3  =  5.0  and  k  =  pi  +  y  =  1.0.  The  integration  time  is  T  =  5  with  an  integration  step  size  of  t  =  0.1  and  a  grid 
resolution  of  0.01  in  both  /  and  p  (momentum).  The  optimal  path  to  extinction  from  the  endemic  state  to  the  disease-free  state 
and  its  counterpart  optimal  path  from  the  disease-free  state  to  the  endemic  state  (found  by  computing  the  backward  FTLE  field) 
are  given  by  the  ridges  of  maximal  FTLE  and  are  overlaid  with  the  analytical  solution  of  the  optimal  paths. 


examples,  T  was  chosen  to  be  sufficiently  large  so  that 
one  obtains  a  measurable  exponential  separation  of  tra¬ 
jectories.  In  figure  3  a,  one  can  see  that  the  optimal  path 
to  extinction  is  given  by  the  ridge  associated  with  the 
maximum  FTLE.  In  fact,  by  overlaying  the  forward 
and  backward  FTLE  fields,  one  can  see  all  three  zero- 
energy  curves  including  the  optimal  path  to  extinction. 
Also  shown  in  figure  3  a  are  the  analytical  solutions  to 
the  three  zero-energy  curves  given  by  p  =  0,  q=  0  and 
equation  (4.5).  There  is  excellent  agreement  between 
the  analytical  solutions  of  all  three  curves  and  the 
ridges  that  are  found  through  numerical  computation 
of  the  FTLE  flow  fields. 

It  is  possible  to  compute  analytically  the  action 
along  the  optimal  path  for  a  range  of  A /pi  values. 
Using  equation  (4.5),  it  is  easy  to  show  that  the 
action  S  is 

5  =2(1  —In  2)-.  (4.6) 

fX 

It  is  clear  from  equation  (4.6)  that  the  action  scales 
linearly  with  A /pi. 

4-2.  Susceptible- infectious -susceptible 
epidemic  model:  external  fluctuations 

We  now  consider  the  well-known  problem  of  extinction 
in  a  susceptible -infectious -susceptible  (SIS)  epidemio¬ 
logical  model,  which  is  a  core  model  for  the  basis  of 
many  recurrent  epidemics.  The  SIS  model  is  described 
by  the  following  system  of  equations: 

S  =  pi  —  piS  +  yl  —  f3IS  (4.7a) 

and 

i=-(fi+y)I  +  piS,  (4.76) 

where  pi  denotes  a  constant  birth  and  death  rate,  (3 
represents  the  contact  rate  and  y  denotes  the  recovery 
rate.  Assuming  the  total  population  size  is  constant 


and  can  be  normalized  to  S+I=  1,  then  equations 
(4.7a,  b)  can  be  rewritten  as  the  following  one¬ 
dimensional  process  for  the  fraction  of  infectious 
individuals  in  the  population: 

i  =  -(n  +  y)I  +  l3I(l-I).  (4.8) 

The  stochastic  version  of  equation  (4.8)  is  given  as 

i  =  -(m  +  y)I  +  j3/(l  -  I)  +  o-ri(t) 

=  F(I)  +  <777  (t),  (4.9) 

where  r\(t)  is  uncorrelated  Gaussian  noise  with  zero 
mean  and  cr  is  the  standard  deviation  of  the  noise  inten¬ 
sity.  The  noise  models  random  migration  to  and  from 
another  population  [15,36]. 

Equation  (4.9)  has  two  equilibrium  points  given  by 
1=  0  (disease-free  state)  and  1=1  —  (pi  +  y)/ (3  (ende¬ 
mic  state).  Using  the  Euler -Lagrange  equation  of 
motion  [42]  from  the  Lagrangian  determined  by 
equation  (4.9)  (£(/,/)  =  [^7(/)]2  =  [/ —  E(/)]2)  along 
with  the  boundary  conditions  given  by  the  extinct 
and  endemic  states,  one  finds  that  the  optimal  path  to 
extinction  (as  well  as  its  counterpart  path  from  the  dis¬ 
ease-free  state  to  the  endemic  state)  is  given  by 
i=±F(i). 

As  in  the  first  example,  one  can  numerically  compute 
the  optimal  path  to  extinction  using  the  FTLE. 
Figure  3  b  shows  the  forward  and  backward  FTLE 
plot  computed  for  T  =  5,  with  (3  =  5.0  and  k  =  pi  + 
y=1.0.  Note  that  we  can  consider  the  combination 
pi  +  y  since  the  Lagrangian  depends  only  on  the  combi¬ 
nation  rather  than  on  either  parameter  separately. 
In  figure  36,  one  can  see  that  the  optimal  path  from 
the  endemic  state  to  the  disease-free  state  is  given  by 
the  ridge  associated  with  the  local  maximum  FTLE. 
Also  shown  in  figure  36  is  the  counterpart  optimal 
path  from  the  disease-free  state  to  the  endemic  state 
(found  by  computing  the  backward  FTLE  field).  In 
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addition,  the  agreement  with  the  analytical  prediction 
is  excellent,  as  shown  in  figure  3b. 

If  one  solves  equation  (4.9)  for  I{t)  and  substitutes 
both  I{t)  and  I{t)  into  the  Lagrangian,  then  one  can 
analytically  find  an  expression  for  the  action  along  the 
optimal  path.  The  expression  for  the  action  is  a  function 
of  k  and  the  reproductive  number  R0  and  is  given  by 


2k(J?0  -  l)3 

3  Rl 


(4.10) 


where  R{)  =  f3 /  k. 


4-3.  Susceptible- inf ectious- susceptible 
epidemic  model:  internal  fluctuations 

We  next  consider  the  one-dimensional  stochastic  ver¬ 
sion  of  the  SIS  epidemic  model  for  a  finite  population 
with  internal  fluctuations  using  the  transition  rates 
found  in  the  electronic  supplementary  material,  §4. 
Using  the  formalism  of  Gang  [35],  one  then  has  the 
following  Hamiltonian  associated  with  this  model: 

H(I,p)  =  (n  +  y)I(e-p - 1)  +  0/(1  - I)(ep - 1),  (4.11) 

where  I  is  the  fraction  of  infectious  individuals  and  p  is 
the  momentum.  Internal  fluctuations  arise  from  the 
random  interactions  of  the  population.  Although  there 
is  no  analytical  solution  for  the  optimal  path  to  extinc¬ 
tion,  we  can  once  again  determine  the  optimal  path  by 
computing  the  FTLE  flow  field  associated  with  this 
system.  Figure  4  shows  the  forward  FTLE  plot 
computed  using  Hamilton’s  equation  of  motion  for 
T  =  10,  with  /3  =  2.0  and  k  =  pi  +  y  =  1.0,  and  as  in 
previous  examples,  the  optimal  path  to  extinction 
from  the  endemic  state  to  the  disease- free  state  is  appar¬ 
ent.  Note  that  the  non-zero  momentum  corresponding 
to  the  extinct  state  qualitatively  agrees  with  similar 
boundary  conditions  found  in  Dykman  et  al.  [9]  and  is 
associated  with  non-zero  probability  flux. 

Once  again,  it  is  possible  to  compute  the  action  along 
the  optimal  path  for  a  range  of  values  of  the  reproductive 
number  R0.  In  contrast  to  the  prior  two  examples,  here 
the  action  must  be  computed  numerically.  Moreover, 
even  the  optimal  path  must  be  found  numerically  using 
the  FTLE  plot  generated  for  each  value  of  R0. 

Given  an  i?0,  we  computed  the  forward  FTLE  flow 
field  using  a  grid  resolution  of  0.005  in  both  position 
and  momentum  space.  To  find  the  optimal  path  corre¬ 
sponding  to  the  ridge  of  maximal  FTLE  values,  we  let 
the  deterministic,  endemic  steady  state  be  the  starting 
point  zq  of  the  path.  Then  a  second  point  z1  on  the 
path  was  found  by  locating  the  maximum  of  the  FTLE 
values  in  an  arc  of  radius  20  grid  points  spanning 
nearly  tt  radians  for  negative  momentum  values  and 
originating  at  zq.  Subsequent  points  zn+1  were  found  by 
taking  the  maximum  FTLE  value  on  an  arc  of  radius 
10  grid  points  spanning  nearly  i t  radians  originating  at 
the  most  recent  point  zn  and  centred  around  the  vector 
zn  —  zn_!  until  the  extinction  line  was  reached.  The  com¬ 
plete  optimal  path  was  estimated  from  the  zn  using  cubic 
spline  interpolation.  Once  the  optimal  path  had  been 
found,  the  action  was  computed  by  numerically  integrat¬ 
ing  the  Lagrangian  along  this  path. 
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Figure  4.  FTLE  flow  field  for  the  SIS  epidemic  model  with 
internal  fluctuations  computed  using  the  Hamiltonian  given 
by  equation  (4.11)  with  /3  =  2.0  and  K=p+  1.0.  The 
integration  time  is  T  —  10  with  an  integration  step  size  of 
t  =  0.1  and  a  grid  resolution  of  0.005  in  both  /  and  p  (momen¬ 
tum).  The  optimal  path  to  extinction  is  given  by  the  ridge  of 
maximal  FTLE. 


Figure  5  a  shows  the  numerically  computed  action 
versus  reproductive  number  over  the  range  1.1  <  R0  <  20. 
The  inset  of  figure  5  a  shows  the  portion  of  figure  5  a 
from  1.1  <  R0  <  2.  Also  shown  in  the  inset  of 
figure  5  a  is  an  analytical,  asymptotic  scaling  result 
that  is  valid  for  values  of  R0  close  to  unity.  As  can  be 
see  in  figure  5  a,  there  is  a  good  agreement  between 
the  numerically  computed  action  and  the  analytically 
computed  action  near  R0=  1.  Details  of  the  derivation 
of  the  analytical  scaling  law  can  be  found  in  the 
electronic  supplementary  material,  §5. 

Figure  5  b  shows  the  numerically  simulated  mean 
extinction  time  versus  reproductive  number  for  the 
one-dimensional  stochastic  SIS  model  for  a  finite  popu¬ 
lation  (see  §4.3).  Also  shown  in  figure  5 b  is  the 
analytical  prediction  found  using  the  previously  men¬ 
tioned  scaling  law  derived  for  values  of  R0  close  to 
unity.  As  one  can  see,  the  agreement  is  excellent. 

5.  CONCLUSIONS 

In  all  of  the  examples  of  §4,  we  have  shown  that  the  opti¬ 
mal  path  to  extinction  is  equated  with  having  a  (locally) 
maximal  sensitivity  to  initial  conditions.  Even  though 
many  possible  paths  to  extinction  exist,  the  dynamical 
systems  approach  converges  to  the  path  that  maximizes 
extinction.  The  parameter  values  chosen  for  the  three 
examples  are  such  that  the  extinct  and  endemic  states 
are  far  away  from  each  other.  Therefore,  in  general, 
there  will  be  no  possible  approximate  analytical  treat¬ 
ment  as  performed  in  Dykman  et  al.  [9].  In  addition, 
we  have  shown  how  to  constructively  compute  numeri¬ 
cally  the  action  for  a  wide  range  of  reproductive 
numbers.  Our  method  allows  for  the  computation  of 
extinction  times  and  can  be  extended  to  high¬ 
dimensional  problems. 

Because  the  method  is  general,  and  unifies  dynami¬ 
cal  systems  theory  with  the  probability  of  extinction, 
we  expect  that  any  system  found  in  other  fields  can  be 
understood  using  this  approach.  Indeed,  in  problems 
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Figure  5.  (a)  Numerically  computed  action  (integrated  along  the  optimal  path  found  numerically  from  the  FTLE  flow  field) 
versus  reproductive  number  R0  for  the  SIS  epidemic  model  with  internal  fluctuations.  The  inset  shows  a  portion  of  the  graph 
near  R0  =  1.  The  numerically  computed  action  is  given  by  the  black  points,  while  the  dashed  curve  shows  an  asymptotic  scaling 
result  that  is  valid  near  R0  =  1,  and  is  given  by  S(R0)  =  (R0  —  1)2/[R0  (1  +  R0)]  (see  the  electronic  supplementary  material, 
§5).  ( b )  Numerically  simulated  (solid  curve  with  black  points)  mean  extinction  time  versus  reproductive  number  for  the  SIS  epi¬ 
demic  model  with  internal  noise  and  the  analytical  prediction  (dashed  curve)  found  using  the  asymptotic  scaling  law  that  is  valid 
near  R0=  1 . 


of  general  extinction,  it  is  now  possible  to  evolve  natu¬ 
rally  to  the  optimal  path,  and  thus  predict  the  path 
that  maximizes  the  probability  to  extinction.  Future 
work  in  this  area  will  include  improved  sampling 
methods  to  find  the  optimal  path  more  efficiently  in 
higher  dimensional  models.  Specific  applications  of 
optimal  path  location  in  the  future  will  include  spatio- 
temporal  extinction  processes  such  as  those  that  occur 
in  pre- vaccine  measles  [7,43]  and  multi-strain  extinction 
in  diseases  such  as  dengue  [44]. 

Finally,  the  optimal  path  method  may  lead  to  novel 
monitoring  and  control  strategies.  In  one  biological 
application,  knowledge  of  the  most  probable  path  to 
extinction  may  help  with  monitoring  an  environment 
and  with  providing  an  estimate  of  the  likelihood  of 
extinction  based  on  where  the  population  lies  relative 
to  the  path.  It  is  known  that  once  a  trajectory  has  left 
the  neighbourhood  of  the  endemic  state,  most  paths 
to  extinction  occur  near  the  optimal  path.  This 
phenomenon  can  be  seen  in  figure  1,  where  the  optimal 
path  lies  on  the  peak  of  the  probability  density  of 
extinction  prehistory.  Therefore,  the  optimal  path  pro¬ 
vides  a  good  location  to  monitor  the  system  for 
possible  extinction  events.  Furthermore,  although  the 
momenta  (noise  effects)  are  not  directly  observable,  it 
may  be  possible  to  infer  them  dynamically  [45]  from 
time  series  data  of  observable  quantities  in  conjunction 
with  the  equations  for  the  time  evolution  of  position 
and  momentum  variables.  Knowledge  of  where  a 
system  lies  in  position-momentum  space  could  provide 
an  estimate  for  how  quickly  a  desired  epidemic  extinc¬ 
tion  could  occur  or  could  provide  the  risk  of  extinction 
for  a  species  one  wishes  to  conserve. 

In  yet  another  application,  knowledge  of  the  optimal 
path  to  extinction  has  the  advantage  that,  in  the  pres¬ 
ence  of  noise  (that  is  estimated  from  data)  and  a  known 
population  of  infectious  individuals,  it  may  be  possible 
to  develop  better  vaccine  controls  that  reduce  the 
time  to  extinction.  Figure  2  a  shows  a  schematic  of  the 
path  to  extinction,  where  the  extinct  state  is  a  saddle 
with  stable  and  unstable  directions.  For  many  epidemic 


models,  the  extinct  state  can  be  shown  to  have  a  similar 
geometry  of  stable  and  unstable  directions.  An  approach 
to  the  extinct  state  on  the  optimal  path  will  lead  to  the 
fastest  time  to  extinction.  Moreover,  since  the  extinct 
state  has  a  saddle  structure  in  the  presence  of  noise,  it 
may  be  controlled  with  projection  methods  [46,47]  or 
probabilistic  techniques  [48]. 

One  can  consider  instituting  a  method  of  parameter 
control,  where  the  parameters  could  be  vaccine  levels  or 
treatment  of  infectious  individuals.  Combined  with  the 
monitoring  techniques  mentioned  above,  the  control 
method  will  allow  one  to  move  an  existing  state  that 
deviates  from  the  optimal  path  closer  to  the  optimal 
path  as  time  evolves.  By  adjusting  the  parameters,  we 
may  target  the  stable  directions  of  the  path  when  we 
are  close  to  epidemic  extinction  [46] .  Comparing  obser¬ 
vations  with  model  predictions  of  the  optimal  path 
allow  us  to  use  controls  to  minimize  the  time  to 
extinction. 
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APPENDIX  A.  THEORY  OF  LARGE 
FLUCTUATIONS 


Letting  A  denote  the  population  size,  the  state  variables 
lGRn  of  the  system  describe  the  components  of  a 
population,  while  the  random  state  transitions,  which 
govern  the  dynamics  are  described  by  the  transition 
rates  kF(X,r),  where  r  E  Rn  is  an  increment  in  the 
change  of  X.  In  the  large  population  limit  without 
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any  fluctuations,  the  mean  field  equations  are  given  by 
the  system  X  =  ^2r  rW(X,r). 

Consideration  of  the  net  change  in  increments  of  the 
state  of  the  system  results  in  the  following  master 
equation  for  the  probability  density  p(X,t)  of  finding 
the  system  in  state  X  at  time  t: 

m^2[W(X-  r;  r)p(X  -  r,  t) 

r 

-W(X;r)p(X,t)\.  (Al) 

The  solution  to  the  master  equation  (A  1)  has  a  peak 
around  a  stable  steady  state  in  the  limit  of  large  popu¬ 
lation  (N^>  1)  with  width  oc  N1^2  [9,23].  However,  since 
we  are  interested  in  the  probability  of  extinction,  we  will 
consider  the  tail  of  the  distribution,  which  gives  the 
probability  of  having  a  small  number  of  individuals. 
The  tail  of  the  distribution  can  be  obtained  by  employ¬ 
ing  the  ansatz  given  by  equation  (2.1),  which  is  an 
assumption  for  how  the  probability  density  scales 
with  population  size  [9,22,49].  Equation  (2.1)  also 
implies  that  a  maximum  in  the  extinction  probability 
can  be  found  minimizing  the  action  over  a  set  of  extinc¬ 
tion  paths  starting  from  the  stable  endemic  state.  The 
assumption  of  this  functional  form  for  p  allows  the 
action  to  be  derived  from  properties  of  the  master 
equation. 

The  density,  p(X,t),  can  be  found  by  substituting 
the  ansatz  given  by  equation  (2.1)  into  the  master 
equation.  The  resulting  equations  for  the  action  are 
given  by  the  Hamilton -Jacobi  equation  for  a 
Hamiltonian,  H ,  given  by  dS/dt  +  H(x,dS/dx;t)  =  0, 
where  we  have  normalized  the  population  components 
and  rates,  respectively,  as  x  =  X/N,  w(x,r)  = 
W(X,r)/N. 

Following  the  classical  mechanics  convention,  define 
a  conjugate  momentum  to  the  state  space,  cc,  by  letting 
p  =  dS/dx  and  where  H(x,p;t)  is  the  classical  Hamil¬ 
tonian  [22].  The  Hamiltonian,  H,  depends  both  on  the 
state  of  the  system,  cc,  as  well  as  the  momentum,  p, 
which  provides  an  effective  force  owing  to  stochastic 
fluctuations  on  the  system.  The  Hamiltonian  equations 
of  motion  provide  trajectories  in  time  of  x{t)  and  p(t), 
and  as  such,  describe  a  set  of  paths  going  from  an  initial 
point  at  time  tt  to  some  final  point  at  time  tf  in  (x,p) 
space.  For  a  given  path,  the  action  is  given  by 
S  =  $lf  p(t)x(t)dt,  and  as  such  determines  the  expo¬ 
nent  of  the  probability  of  observing  that  path.  (It 
should  be  noted  that  instead  of  using  the  Hamiltonian 
representation,  one  could  use  the  Lagrangian  represen¬ 
tation  T(cc,  x\  t )  =  —H{x,p\t)-\-x-  p ,  which  results  in 
an  equivalent  solution.) 

The  action  reveals  much  information  about  the  prob¬ 
ability  evolution  of  the  system,  from  scaling  near 
bifurcation  points  in  non- Gaussian  processes  to  rates 
of  extinction  as  a  function  of  epidemiological  par¬ 
ameters  [9,49].  As  already  stated,  in  order  to 
maximize  the  probability  of  extinction,  one  must  mini¬ 
mize  the  action  S.  The  minimizing  formulation  entails 
finding  the  solution  to  the  Hamilton- Jacobi  equation, 
which  means  that  one  must  solve  the  2  n- dimensional 
system  of  Hamilton’s  equations  [x  =  dpH(x,  p), 


p  =  —dxH(x.  p)]  for  x  and  p,  where  the  Hamiltonian 
is  explicitly  given  as 

H(x,p,  t)  =  ^2  w(x,  r)[exp(pr)  -  1].  (A  2) 

r 

The  appropriate  boundary  conditions  of  the  system 
are  such  that  a  solution  starts  at  a  non-zero  state, 
such  as  an  endemic  state,  and  asymptotically 
approaches  one  or  more  zero  components  of  the  state 
vector,  representing  a  disease-free  state.  Therefore,  a 
trajectory  that  is  a  solution  to  the  two-point  boundary 
value  problem  determines  a  path,  which  in  turn  yields 
the  probability  of  going  from  the  initial  state  to  the 
final  state.  The  optimal  path  to  extinction  is  the  path 
that  minimizes  the  action  in  either  the  Hamiltonian  or 
Lagrangian  representation. 

We  compute  the  trajectory  satisfying  the  Hamil¬ 
tonian  system  that  has  as  its  asymptotic  limits  in 
time  the  endemic  state  as  t^—oo  and  the  extinct 
state  as  t  ^  +  oo.  The  momentum  p  represents  the 
force  of  the  fluctuations  on  the  population,  and  this 
momentum  changes  the  stability  of  the  equilibrium 
points.  Both  the  endemic  and  extinct  states  have 
attracting  and  repelling  directions  for  p  ^  0,  as  shown 
schematically  in  figure  2. 

Given  optimal  path  trajectories  (&(£),  p(t)),  the  action 
with  correct  limits  is  found  by  S (x)  =  pxdt  [42]. 
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